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ABSTRACT 

A  model  of  formation,  development,  and  acoustic  properties  of  the  bubble  cloud  resulting  from  an  underwater  ex¬ 
plosion  is  presented.  The  model  includes  several  parts:  explosion  globe  dynamics,  initial  break-up  of  the  explosion 
globe,  turbulence  created  by  the  explosion  globe  fragmentation,  and  further  break-up  of  the  bubbles  by  the  turbu¬ 
lence.  The  time  history  of  the  bubble  cloud  properties  is  calculated  under  the  assumption  of  the  cloud  being  a  collec¬ 
tion  of  non-interacting  bubbles.  Model  results  are  compared  with  the  available  experimental  data. 


INTRODUCTION 

It  is  of  interest  for  some  defence  applications  to  understand 
the  dynamics  and  to  model  the  acoustic  properties  of  the 
bubble  cloud  remnant  of  the  underwater  explosion.  There  are 
not  many  papers  published  on  the  subject.  We  are  only  aware 
of  one  (Holt&Culver,  2011),  in  which  the  remnant  bubble 
cloud  is  investigated  experimentally  by  acoustic  method.  To 
the  best  of  our  knowledge,  there  were  no  attempts  to  develop 
a  theoretical  model  of  the  remnant  bubble  cloud  formation 
and  dynamics.  The  present  document  reports  the  develop¬ 
ment  of  such  a  model.  We  apply  a  physics-based  approach 
and  avoid  complex  and  computationally  demanding  numeri¬ 
cal  simulations.  Such  an  approach  will  undoubtedly  require 
certain  approximations  and  the  introduction  of  some  empiri¬ 
cal  parameters  into  the  model.  These  will  have  to  be  adjusted 
from  the  comparison  with  the  existing  or  future  experimental 
data  or  high-fidelity  numerical  simulations.  We,  however,  are 
not  aware  of  any  such  accurate  numerical  simulations  of  this 
problem.  Taking  into  account  the  complexity  of  the  phe¬ 
nomenon,  it  is  unlikely  that  such  simulations  could  be  per¬ 
formed  in  the  near  future. 

MODEL  OUTLINE 

As  a  result  of  an  underwater  explosion  (UNDEX),  a  gas 
globe  is  formed,  which  experiences  several  expansions  and 
contractions  during  short  time  after  the  detonation.  We  will 
use  the  term  ’’globe”  for  the  initial  gas  bubble  formed  by  the 
underwater  explosion  to  distinguish  it  from  the  bubbles  as 
product  of  the  initial  globe  disintegration.  During  these  oscil¬ 
lations  the  explosion  globe  rises  in  water,  especially  quickly 
during  contraction  phases  when  the  drag  force  is  lower.  After 
two  to  three  oscillations  the  large  explosion  bubble  disinte¬ 
grates  into  a  large  number  of  smaller  bubbles,  which  then  rise 
towards  the  surface  of  water  with  various  speeds  depending 
on  their  size.  The  purpose  of  this  model  is  to  estimate  the 
bubble  size  distribution,  the  bubble  spatial  distribution  and 
the  dynamics  of  their  rise.  The  acoustic  properties  of  the 
remnant  bubble  cloud  can  be  easily  derived  from  the  known 
bubble  size  distribution. 

The  steps  of  the  currently  suggested  model  are  summarised 
below.  The  references  to  specific  model  elements  will  be 
given  in  the  corresponding  sections  of  the  paper. 


-  The  oscillations  and  rise  of  the  explosion  globe  are  based  on 
the  models  available  in  the  literature,  with  slight  modifica¬ 
tions,  of  spherical  bubble  shape  dynamics  and  motion.  There 
are  more  sophisticated  numerical  models  describing  the 
UNDEX  globe  oscillation  and  rise,  but  for  the  purpose  of  this 
research  the  current,  computationally  efficient  approach  will 
suffice. 

-  It  is  assumed  that  at  the  third  minimum  of  the  UNDEX 
globe  oscillation,  it  disintegrates  into  smaller  bubbles.  The 
size  distribution  of  the  fragments  is  estimated  from  the 
growth  rate  of  the  modes  of  the  Rayleigh-Taylor  instability 
for  spherical  cavity.  It  is  assumed  that  initial  pressure  of  gas 
in  the  fragments  is  the  same  as  in  the  explosion  globe  just 
before  the  break-up.  Then  the  bubbles  expand  to  bring  the 
internal  gas  pressure  into  the  balance  with  the  ambient  pres¬ 
sure.  The  resulting  bubble  cloud  size  is  estimated  as  being 
proportional  to  the  total  gas  volume.  The  coefficient  of  pro¬ 
portionality  is  treated  as  an  empirical  parameter  in  the  current 
model. 

-  The  explosion  globe  potential  energy  at  the  minimum  of  its 
radius  goes  into  the  kinetic  energy  of  the  turbulent  spot  in  the 
fluid.  The  time-space  distribution  of  the  turbulent  kinetic 
energy  is  obtained  by  solving  the  corresponding  equation  in 
approximation  of  spherical  symmetry. 

-  The  bubbles  resulting  from  the  initial  fragmentation  of  the 
explosion  globe  are  further  broken-up  by  turbulence.  To  de¬ 
scribe  this,  one  of  the  bubble  break-up  models  is  used.  It 
should  be  noted  here  that  all  of  the  bubble  break-up  models 
by  turbulence  are  developed  in  the  assumption  of  small  gas 
volume  fraction  and  isotropic  nature  of  turbulence.  In  this 
problem,  the  gas  volume  fraction  is  high  and  the  turbulence  is 
not  isotropic.  However,  in  the  absence  of  a  better  model,  one 
of  the  existing  models  is  applied  with  a  note  of  addressing 
this  problem  in  future.  The  unsteady  equation  of  the  bubble 
population  balance  is  solved  using  the  time-space  distribution 
of  the  turbulence  obtained  in  the  previous  step. 

-  The  rise  of  the  bubbles  to  the  surface  is  then  calculated  from 
the  balance  of  the  buoyancy  and  drag  forces.  The  spatial  and 
size  distribution  of  the  gas  volume  fraction  is  then  used  to 
estimate  the  acoustical  properties  of  the  gas  bubble  cloud. 
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OSCILLATIONS  AND  RISE  OF  THE 
EXPLOSION  BUBBLE 

There  exist  sophisticated  numerical  models  of  the  explosion 
globe  rise  and  oscillations.  Here,  for  the  sake  of  UNDEX 
remnant  bubble  cloud  model  development,  we  need  a  rela¬ 
tively  simple,  computationally  efficient  estimation  of  the 
bubble  radius  and  pressure  inside  the  bubble  at  the  beginning 
of  the  explosion  bubble  break-up.  In  the  current  model  we 
base  our  estimation  on  a  relatively  simple  equation  of  a 
spherical  bubble  motion  and  oscillations.  Later,  these  estima¬ 
tions  may  be  obtained  from  a  more  sophisticated  numerical 
model. 

The  models  of  the  oscillations  and  motion  of  the  explosion 
bubble,  sometimes  referred  to  as  globe  here  to  distinguish 
from  the  bubbles  resulting  from  the  explosion  globe  fragmen¬ 
tation,  date  back  to  1940s  {Underwater  Explosion  Research, 
1950;  Cole,  1948).  The  early  models  were  based  on  the  as¬ 
sumption  of  incompressible  flow  and  did  not  take  into  ac¬ 
count  the  energy  loss  due  to  pressure  wave  radiation  during 
the  collapse-rebound  stage  of  bubble  oscillation.  As  a  result 
they  did  not  account  for  the  damping  of  bubble  shape  oscilla¬ 
tions  observed  in  experiments.  Recently  a  model  has  been 
published  (Geers&Hunter,  2002),  which  accounts  for  the 
wave  effect.  The  model  is  based  on  the  doubly  asymptotic 
approximation  and  combines  relative  simplicity  based  on  the 
assumption  of  spherical  shape  of  the  explosion  globe  with  the 
reasonable  prediction  of  the  globe  oscillation  damping.  In  our 
model  of  formation  of  the  UNDEX  remnant  bubble  cloud  we 
employ  this  model  of  the  explosion  bubble  oscillations  with 
some  modifications.  Thus  the  radius  of  the  spherical  explo¬ 
sion  globe,  R  ,  is  described  by  the  following  equation: 

RRfl+l^f2+Rf^EfA  +  \^L  (1) 

2  pj  4 

where  pl  is  the  density  of  gas  inside  the  bubble  and 
surrounding  liquid,  respectively,  cg ,  ct  is  the  speed  of  sound 
in  gas  and  liquid,  p ^  is  the  ambient  pressure  of  surrounding 
fluid,  p  is  the  initial  gas  pressure  inside  the  bubble.  We  do 
not  reproduce  here  the  functions  fl9  f2,f3  and  f4  referring 

to  Geers&Hunter  (2002)  (Equations  32  and  33).  The  second 
term  on  the  right-hand  side  of  the  equation  (1)  is  due  to  the 
added  pressure  on  the  bubble  surface  resulting  from  the  bub¬ 
ble  motion  (Hsiao  et  al.,  2003).  We  include  this  term  in  addi¬ 
tion  to  the  model  provided  by  Geers&Hunter  (2002). 

Initial  conditions  for  equation  (1)  are  obtained  from  combina¬ 
tion  of  similitude  relation  for  the  far-field  shock-wave  pres¬ 
sure  profiles  and  the  volume  acceleration  model 
(Geers&Hunter,  2002).  The  bubble  translation  model  de¬ 
scribed  in  the  same  paper  does  not,  however,  provide  a  satis¬ 
factory  result  unless  an  artificial  drag  coefficient  is  intro¬ 
duced.  Therefore,  we  have  chosen  to  use  the  following  equa¬ 
tion  often  applied  to  modelling  of  bubble  motion 
(Hsiao&Pauley,  1999): 


In  this  equation  is  the  velocity  of  the  bubble,  u  is  the 
velocity  of  the  surrounding  fluid,  g  is  the  acceleration  due  to 
gravity,  and  CD  is  the  drag  coefficient.  Here  we  use  the 
Grace  Drag  model  (Clift  et  al.,  1978;  ANSYS  CFX-Solver , 
2010),  which  accounts  for  wide  range  of  bubble  Reynolds 
number  and  various  regimes  of  bubble  deformation. 


Figure  1.  Explosion  globe  radius  (top  plot)  and  depth  (bot¬ 
tom  plot). 

Here  we  will  consider  an  example  from  the  laboratory 
UNDEX  experiment  conducted  at  the  Carderock  Division  of 
the  Naval  Surface  Warfare  Center  (Harris  et  al.,  2010).  In  this 
experiment  0.23  kg  (0.5  lb)  of  pentolite  explosive  was  deto¬ 
nated  at  2.29  m  (7.5  feet)  of  depth.  The  result  of  applying  the 
above  model  to  this  case  is  presented  in  Figure  1  together 
with  the  experimental  results  obtained  by  processing  video 
taken  during  the  experiment.  It  should  be  noted  here  that  the 
above  mentioned  report  does  not  specify  the  exact  composi¬ 
tion  of  explosive  used  in  the  experiment  and  there  are  several 
versions  of  the  explosive  called  pentolite.  In  the  model  reali¬ 
sation  we  used  the  physical  parameters  for  pentolite  available 
in  the  literature.  More  accurate  data  on  the  physical  proper¬ 
ties  of  the  specific  explosive  used  in  the  test  may  improve  the 
agreement  between  the  model  and  experimental  data.  Never¬ 
theless,  the  current  model  gives  much  better  agreement  than 
the  usual  bubble  oscillation  models  which  do  not  take  into 
account  the  energy  loss  due  to  shock  wave. 


3  1/  \i  |  3  /  \dR 
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EXPLOSION  GLOBE  FRAGMENTATION 

We  assume  that  the  explosion  globe  break-up  occurs  at  the 
third  minimum  as  a  result  of  instability  of  the  bubble  spheri¬ 
cal  shape.  It  is  a  simplification  because  in  reality  the  frag¬ 
mentation  of  the  explosion  globe  occurs  in  steps  at  the  sec¬ 
ond,  third,  and  even  fourth  minimum  of  the  breathing  modes, 
which  could  be  seen  from  the  video  sequence  (Harris  et  al., 
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2010).  However  we  believe  that  compressing  this  process 
into  one  single  act  of  fragmentation  is  a  reasonable  assump¬ 
tion  which  significantly  simplifies  the  model.  The  model 
improvement  in  this  respect  may  be  addressed  in  future  work. 

Just  before  the  explosion  globe  break-up  its  radius  is  R*  and 
the  gas  pressure  inside  the  bubble  .  In  the  current  model 

we  obtain  these  values  from  the  solution  of  equations  in  the 
previous  section.  Alternatively,  they  can  be  obtained  from  a 
more  accurate  numerical  model  of  underwater  explosion. 

We  estimate  initial  bubble  size  distribution  from  the  modes 
of  the  Rayleigh-Taylor  instability  for  the  spherical  bubble. 
The  growth  rate  of  the  non-spherical  distortions  on  the  bub¬ 
ble  surface  is  proportional  to  the  following  value  (Brennen, 
2002): 

f(m)  =  (m  -  l)[r  -  (m  +  \\m  +  2)]  (3) 

where  r  =  pR2R/  cr  and  m  is  the  order  of  a  spherical  har¬ 
monic  distortion.  We  postulate  that  each  mode  m  leads  to 
bubbles  of  size  rm  =  R*0  / m  in  the  daughter  bubble  cloud 

and  that  the  number  of  bubbles  of  this  size  is  proportional  to 
the  function  f(m)- 

nm  =  A(m  -  l)[r*  ~(m  + 1  \m  +  2)] 


f  Rl 

f  r!  ] 

K 

-1 

r0*  - 

-5-+1 

— +  2 

Vr» 

) 

) 

where  ro*  is  the  value  of  T  at  R  =  R*.  The  normalisation 

coefficient  A  is  determined  from  the  balance  of  the  gas  vol¬ 
ume: 

(5) 


We  assume  that  straight  after  the  break-up  the  gas  in  the 
daughter  bubbles  has  the  same  pressure,  p *Q ,  as  the  explo¬ 
sion  bubble  just  before  break-up.  Then  the  bubbles  expand  to 
reach  the  balance  with  the  pressure  in  the  surrounding  fluid. 
We  assume  here  that  during  globe  break-up  the  gas  tempera¬ 
ture  rapidly  adjusts  to  the  temperature  of  water  and  subse¬ 
quent  expansion  of  break-up  products  is  isothermal.  The  new 
bubble  size  can  be  obtained  from  the  following  equation: 

P^=p2V 2,  (6) 

where  px  =  p*g0  and  p2  =  pa  +  pghb +2cr /rm,  hb  being 
the  depth  of  the  bubble,  and  pa  the  atmospheric  pressure. 

This  results  in  the  following  cubic  equation  for  the  new  bub¬ 
ble  radius,  r  : 


r^3  2<J  2 

r  H - r  —  - 


Pg  o 


Pa  +  Pghb 


Pa  +  Pg\ 


-r:  =  O' 


(7) 


If  the  surface  tension  can  be  neglected,  this  equation  simpli¬ 
fies  significantly: 


however,  with  the  modem  computers  solving  the  cubic  equa¬ 
tion  does  not  affect  the  overall  computation  time  noticeably. 
The  new  bubble  size  distribution  is  calculated  as: 


(9) 


and  the  new  total  volume  of  the  bubbles  in  the  cloud  is: 


The  radius  of  the  bubble  cloud  can  be  estimated  as  follows: 


R-af 


/3^V/3 


An 


(11) 


Currently  we  consider  ac  as  an  empirical  parameter.  Obvi¬ 
ously,  ac>  1  • 


Continuing  with  the  example  from  the  previous  section  and 
assuming  ac- 1.5,  we  obtain  the  bubble  size  distribution 

shown  in  Figure  2. 


Figure  2.  Estimated  bubble  size  distribution  after  the  initial 
break-up  of  the  explosion  bubble. 


TURBULENCE  RESULTING  FROM  EXPLOSION 
BUBBLE  BREAK-UP 

We  assume  that  almost  all  potential  energy  of  the  explosion 
bubble  before  break-up  transfers  into  the  turbulent  kinetic 
energy,  which  then  gradually  decays  due  to  dissipation, 
breaking  further  the  bubbles  in  the  remnant  bubble  cloud. 
The  potential  energy  of  the  compressed  explosion  globe  is 
computed  as  the  work  performed  against  external  pressure 
during  globe  adiabatic  expansion  from  the  compressed  state 
to  the  state  of  equilibrium  with  the  external  pressure  (Cole, 
1948): 


i 

IT  =  4zr  J  [ \p{r)~  P^\2dr' 


(12) 


r  *  y/3 

Pg  o  , 

Pa+PSK, 


Here  =  pa  +  pghb  is  the  ambient  pressure  at  the  bubble 
depth,  hb  •  During  adiabatic  expansion  the  pressure  inside  the 
bubble  is  changing  according  to  the  following  equation: 
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(13) 


Obviously,  R  is  obtained  from  equation  (13)  when 
p(R{)  =  Prjo  .  Integration  of  (12)  will  give  the  following  equa¬ 
tion  for  the  bubble  potential  energy: 


4  71 


i -y 


f  P  A 

fo 


-1 


-P 


-1 


(14) 


When  bubble  breaks  its  potential  energy  mainly  goes  into  the 
kinetic  energy  of  the  fluid  creating  a  turbulent  spot  with  the 
total  kinetic  energy 

*0  =  An>  (15) 

where  j3K  <  1  accounts  for  energy  loss  due  to  change  in  sur¬ 
face  energy  and  acoustic  radiation.  In  the  current  research  we 
neglect  these  losses  and  assume  j3k  «  1 .  We  assume  also  that 
the  initial  radius  of  the  turbulent  spot  coincides  with  the  ra¬ 
dius  of  the  bubble  cloud,  R  .  Thus  the  initial  condition  for  the 
turbulent  kinetic  energy  per  unit  mass  is: 

*  =  {*»>  (16) 
[o,  r  >  R 


where  ^  _  ^K0  .  The  decay  of  the  turbulent 

0  4^(1  -a-3)R3 

kinetic  energy  is  described  by  the  following  equation  (Wil¬ 
cox,  1994): 


dJL__c  fc3/2  f  1  8  f Jkl^dk 

dt~  M  l  +  r2  dr  a,  r  dr 


(17) 


For  the  turbulent  parameters  in  this  equation  we  use  the  fol¬ 
lowing  standard  values:  CM  =0.09,  <Jk  =  \,  l  =  2R  ■  Equation 

(17)  is  solved  numerically  using  the  Crank-Nicolson  finite 
difference  method.  The  dissipation  rate  of  the  turbulent  ki¬ 
netic  energy,  which  is  used  further  in  the  model  of  the  bubble 
break-up,  is  calculated  according  to  the  following  equation: 

r  3/2 

£  =  CU—  08) 

Continuing  with  the  example,  Figure  3  shows  the  decay  of 
the  turbulent  kinetic  energy  at  the  centre  of  the  turbulent  spot 
obtained  from  the  numerical  solution  of  equation  (17). 


Figure  3.  Decay  of  the  turbulent  kinetic  energy  at  the  centre 
of  the  turbulent  spot. 


FURTHER  BUBBLE  BREAK-UP  BY 
TURBULENCE 

The  change  of  the  bubble  size  distribution,  n  ,  due  to  bubble 
break-up  by  turbulence  is  described  by  the  following  equa¬ 
tion  (Lasheras  et  al.,  2002): 


=  ]m(D0)f(D,D0)g(DMD0)dD0  - g{D)n  •  (19) 

ot  D 


The  first  term  on  the  right-hand  side  of  this  equation  de¬ 
scribes  the  birth  rate  of  bubbles  of  diameter  D ;  the  second 
term  is  the  death  rate  of  bubbles  of  this  size.  In  this  equation, 
g(Z))  is  the  break-up  frequency,  /(/),D0)  is  the  probability 

density  function  of  the  daughter  bubbles  resulting  from  the 
break-up  of  the  mother  bubble  of  size  D0,  m(D0 )  is  the  av¬ 
erage  number  of  daughter  bubbles  per  one  break-up  event. 
Here  we  use  the  break-up  frequency  and  daughter  bubbles 
probability  density  function  from  (Martinez-Bazan  et  al., 
1999a,  1999b): 

g(s, D0 )  =  K,D-'  J/?UO0)2  3-l2-^-  ,  (20) 


fip,  D0)=  B 


iPfi(sDr 


Dn 


i  — 


Dn 


(21) 


In  the  above  equations  s  is  the  dissipation  rate  of  the  turbu¬ 
lent  kinetic  energy,  Kg  and  ft  are  empirical  constants, 

D2  =  (l  -  D3J  3  is  the  size  of  the  second  daughter  bubble, 

and  B  is  the  normalisation  factor.  The  following  values  for 
the  empirical  constants,  K  «  0.25  and  fj  =  8.2,  are  sug¬ 
gested  by  Martinez-Bazan  et  al.  (1999a,  1999b). 

Applying  the  turbulent  kinetic  energy  obtained  in  the  previ¬ 
ous  section  to  the  bubble  cloud  with  the  initial  bubble  size 
distribution  from  Figure  2  we  obtain  the  final  bubble  size 
distribution  in  the  bubble  cloud  shown  in  Figure  4. 
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Figure  4.  Final  bubble  size  distribution  after  break-up  by 
turbulence. 

RISE  OF  THE  BUBBLE  CLOUD 

Dissipation  of  the  turbulent  kinetic  energy  in  the  turbulent 
spot  caused  by  initial  break-up  of  the  explosion  bubble  oc¬ 
curs  relatively  quickly.  Therefore,  we  assume  that  in  this  time 
small  bubbles  in  the  bubble  cloud  do  not  move  significantly 
from  their  initial  positions.  Therefore,  we  assume  that  the  rise 
of  the  bubble  cloud  starts  after  the  bubble  break-up  by  turbu¬ 
lence  is  finished.  Such  an  approximation  simplifies  the  prob¬ 
lem  considerably  without,  we  believe,  significant  loss  of 
accuracy.  Another  simplifying  assumption  we  make  is  that 
the  bubbles  rise  at  their  terminal  velocity,  U ,  obtained  from 
the  balance  of  buoyancy  and  drag  force: 

(pl-p^bg~\PlCDU1Ab  (22) 

Here  pb  is  the  density  of  the  gas  inside  the  bubble, 
y  =—7rR3  is  the  bubble  volume,  CD  is  the  drag  coefficient, 

b  3 

and  Ab  =  7rR2  is  the  bubble  projected  area.  Here  again  we 
use  the  Grace  Drag  model  described  above. 

For  each  size  in  the  bubble  size  distribution  we  calculate  the 
bubble  depth,  hb(t),  as  a  function  of  time  from  the  following 

kinematic  equation: 


a  relevant  position  of  their  starting  point  in  the  bubble  cloud 
for  subsequent  computation  of  the  volume  fraction  in  the 
rising  bubble  cloud.  In  this  assumption  equation  (24)  can  be 
written  in  the  following  form: 

(w6>ZZ,)  =  (Wo>Z0  +/z(wi0)^))’  (25) 

where  zc  is  the  z-coordinate  of  the  initial  bubble  cloud  cen¬ 
tre. 


at 


Figure  5.  Logarithm  of  volume  fraction  at  different  time 
(23)  moments  (from  top  down):  1,  5,  and  20  sec  from  start  of  bub¬ 

ble  cloud  rise. 


The  equation  (22)  and  (23)  are  solved  numerically.  In  this 
solution  we  take  into  account  the  change  in  the  bubble  vol¬ 
ume,  and  as  a  result  all  other  bubble  parameters,  due  to 
changing  pressure  of  the  ambient  fluid. 

Under  the  above  assumptions  all  bubble  trajectories  will  be 
vertical  lines  with  the  different  variation  of  bubble  depth  with 
time  depending  on  the  initial  bubble  size.  Mathematically,  we 
can  describe  the  bubble  trajectories  by  the  following  equa¬ 
tions  in  the  coordinate  system  ( x,y,z )  with  x,y  being  hori¬ 
zontal  coordinates  and  z  vertical: 

(xb  ,yb’zb)=(xo’y0’  4zo  >  rm\4)  •  (24) 

While  computing  bubble  trajectories  for  different  bubble 
sizes  we  assume  that  they  all  originate  from  the  centre  of  the 
bubble  cloud.  The  trajectories  are  then  simply  translated  into 


To  calculate  the  volume  fraction  in  the  rising  bubble  cloud 
we  divide  the  initial  bubble  cloud  into  many  volume  elements 
and  originate  bubble  trajectories  for  all  bubble  sizes  from 
each  volume  element.  Each  bubble  trajectory  thus  associated 
with  a  certain  gas  volume.  Obviously,  for  smooth  solution  a 
sufficient  number  of  volume  elements  are  required.  In  this 
example  we  use  100x100x100  volume  elements  for  initial 
volume  subdivision.  Also  a  sufficient  number  of  bubble  size 
groups  are  neded  to  obtain  a  smooth  solution  along  the  verti¬ 
cal  coordinate.  In  this  example  20  size  groups  are  used.  The 
spatial  region  above  the  initial  bubble  cloud  is  divided  into 
volume  elements  as  well  and  from  analysis  of  interjections  of 
bubble  trajectories  with  these  volume  elements  the  volume 
fraction  of  the  bubble  cloud  at  a  certain  time  is  calculated. 
Figure  5  shows  an  example  of  the  volume  fraction  in  the 
rising  bubble  cloud  at  the  different  time  moments. 
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Here  we  should  note  that  for  calculation  of  the  volume  frac¬ 
tion  and  acoustical  properties  of  the  bubble  cloud  we  do  not 
use  the  whole  final  bubble  size  distribution  (Figure  4).  Some 
of  the  bubbles  in  this  distribution  are  still  large  and  will  rise 
to  the  surface  very  quickly.  In  this  research  we  are  interested 
in  the  properties  of  the  bubble  cloud  at  the  relatively  long 
time  after  the  detonation,  on  the  order  of  several  minutes. 
This  means  a  bubble  cloud  consisting  of  small  bubbles. 
Therefore,  in  our  model  we  only  consider  the  part  of  the  final 
bubble  size  distribution  with  bubbles  less  than  a  certain  size. 
In  the  example  here  we  used  the  bubble  diameter  of  5  mm  as 
the  cut-off  bubble  size. 


Once  the  volume  fraction  of  the  bubble  cloud  is  known  for 
the  different  size  bins,  the  backscattering  cross  section  per 
unit  volume  Shs  and  the  spatial  attenuation  rate  ab  due  to 

bubbles  are  given  by  the  following  equations  (Med- 
win&Clay,  1998): 

v  =  Zw<AcrT’/)’  (25) 

1=4 

ai=4.34S4=4.34f>,(rl,/H>  (26) 


where  S  is  the  extinction  cross  section  per  unit  volume.  The 
differential  scattering  cross  section  of  a  single  bubble 


>/)= 


[(/*//)2-l]2+<52 


(27) 


and  the  extinction  cross  section  of  a  single  bubble 


A>/)= 


4;ar2  (<?  /  ) 

lfRlf)2-\]+52 


where  the  total  damping  constant 


(28) 


Figure  6.  Backscattering  cross-section  per  unit  volume  at 
different  time  moments  (from  top  down):  1,  5,  and  20  sec 
from  start  of  bubble  cloud  rise. 


8  =  8r+8t+8v 


(29) 


The  backscattering  cross-section  per  unit  volume  is  computed 
from  the  raw  sonar  data,  S ,  as  follows: 


is  the  sum  of  the  reradiation  term  Sr ,  the  thermal  damping 
term  8t ,  and  the  viscous  damping  term  Sv  •  The  damping 
terms  and  the  resonance  frequency,  fR ,  depend  on  the  physi¬ 
cal  properties  of  the  bubble  gas  and  the  ambient  fluid  as  well 
as  bubble  depth.  The  corresponding  equations  are  given  in 
section  8.2  of  Medwin&Clay  (1998)  and  we  do  not  reproduce 
them  here. 

Example  of  the  backscattering  cross-section  per  unit  volume 
is  shown  in  Figure  6. 

To  compare  the  model  with  available  experimental  data  (Har¬ 
ris  et  al.,  2010),  we  compute  the  mean  backscattering  cross- 
section  per  unit  volume  in  a  horizontal  slice  located  at  the 
depth  of  0.9  m  (three  feet),  or  1.37  m  (4.5  feet)  above  the 
point  of  detonation.  In  the  experiment  the  acoustic  measure¬ 
ments  were  taken  at  this  horizontal  slice  with  Reson  7125 
multibeam  sonar.  Example  of  sonar  raw  data  is  shown  in 
Figure  7.  The  light  coloured  spot  in  the  sonar  data  indicates 
the  cross-section  of  the  bubble  cloud  rising  from  the  point  of 
detonation. 


Sv  =  20 log10  S —  G +  C  —  SL  +  TL  —  N) log10  V  >  (30) 

where  SL  is  the  sonar  source  level  in  dB, 
TL  =  2 aR  +  401og10  R  is  the  two-way  transmission  loss,  C 
is  the  sonar  calibration  factor,  V  is  the  insonified  volume, 
and  G  =  Gv  +  Gc  =  2a0R  +  2/?log10  R  +  Gc  is  the  sonar 

gain,  consisting  of  variable  and  constant  parts.  From  sonar 
data  taken  in  the  air  it  can  be  concluded  that  with  high  degree 
of  accuracy  p  «  20  •  Unfortunately,  the  sonar  was  not  cali¬ 
brated  for  this  experiment  and  the  calibration  constant,  C ,  is 
unknown.  Here  we  attempted  to  calibrate  sonar  using  data  on 
the  gas  volume  fraction  and  bubble  size  distribution  for  the 
bubble  maker,  which  was  also  used  in  the  UNDEX  experi¬ 
ment  (Harris  et  al.,  2010).  The  sonar  data  from  the  bubble 
maker  jet  was  obtained  by  the  same  sonar.  On  the  other  hand, 
the  data  on  the  gas  volume  fraction  and  the  bubble  size  dis¬ 
tribution  in  the  bubble  maker  jet  was  acquired  in  a  different 
laboratory  experiment  using  different  measurement  tech¬ 
niques  by  Dynaflow,  Inc.  Comparing  the  volume  fraction 
data  obtained  from  the  sonar  data  with  the  data  provided  by 
Dynaflow  we  calculated  the  calibration  factor.  The  calcula¬ 
tion  of  the  calibration  factor  was  performed  in  five  different 
locations  of  the  bubble  maker  jet.  The  difference  between 
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maximum  and  minimum  values  of  the  calibration  factor  was 
18  dB.  There  are  many  factors  which  can  contribute  to  the 
uncertainty  in  the  value  of  the  calibration  factor  with  such  a 
non-traditional  approach  to  calibration.  The  most  important 
among  them  is  turbulent  fluctuations  of  the  bubble  maker  jet. 
Other  factors  include  inaccuracy  in  the  measurement  of  the 
volume  fraction  and  the  bubble  size  distribution. 


0.5  1  1.5  2  2.5 


Figure  8.  Backscattering  cross-section  per  unit  volume  calcu¬ 
lated  from  sonar  data  in  a  horizontal  slice  through  the  bubble 
cloud. 

Using  thus  obtained  calibration  factor  we  processed  a  se¬ 
quence  of  sonar  pings  in  a  horizontal  slice  through  the  bubble 
cloud.  Example  of  the  backscattering  cross-section  obtained 
from  the  processing  of  the  sonar  data,  zoomed  onto  the  bub¬ 
ble  cloud  is  shown  in  Figure  8. 

The  black  curve  in  Figure  9  shows  the  mean  backscattering 
cross-section  in  the  horizontal  slice  of  the  bubble  cloud  ob¬ 
tained  from  processing  800  pings.  The  error  bars  on  this 
curve  indicate  the  uncertainty  in  the  calibration  factor.  The 
red  curve  in  this  figure  is  obtained  from  the  model. 


Figure  9.  Backscattering  cross-section  at  the  horizontal  slice 
of  the  rising  bubble  cloud.  Depth  of  the  slice  location  is  3 
feet. 


CONCLUSION 

In  this  paper  we  have  described  the  development  of  the 
model  of  formation,  dynamics  and  acoustic  properties  of  the 
bubble  cloud  resulting  from  an  underwater  explosion.  The 
overall  model  consists  of  several  sub-models  for  different 
steps  of  the  bubble  cloud  formation  and  dynamics.  These 
steps  include  initial  explosion  globe  oscillation  and  rise,  its 
fragmentation  into  the  smaller  bubbles  and  further  break-up 
of  latter  by  the  turbulence  created  by  the  explosion  globe 
break-up.  The  time  history  of  the  bubble  cloud  volume  frac¬ 
tion  and  its  acoustic  properties  is  calculated  under  the  as¬ 
sumption  that  the  bubbles  in  the  cloud  are  not  interacting 
with  each  other  neither  hydrodynamically  nor  acoustically. 
Comparison  of  the  model  with  a  laboratory  experiment  shows 
a  reasonable  agreement.  However,  the  value  of  such  compari¬ 
son  is  somewhat  diminished  by  the  fact  that  the  sonar  used  in 
the  laboratory  measurement  was  not  properly  calibrated, 
which  introduced  a  large  uncertainty  in  the  experimental 
data.  Further  model  improvement  and  validation  will  be  ad¬ 
dressed  in  future. 
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